Hidden multiferroic order in graphene zigzag ribbons 
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The insulating magnetic phase in graphene zigzag ribbons, predicted both by density functional 
and mean field Hubbard model calculations, is described without additional approximations with 
a BCS wave function of two phase-locked condensates of spin-polarized electron-hole pairs. The 
associated order parameter is the spin dipole operator that features both magnetic and electric 
order and accounts for the spin-resolved ferroelectricity of the system. Each condensate is associated 
to a spin-dependent dipole and their relative phase locking sets the total electric dipole and total 
magnetization equal to zero. 

PACS numbers: 



The quest of novel electronic phases, characterized by 
new order parameters, and the interplay between elec- 
tric and magnetic degrees of freedom are two of the 
major themes in condensed matter physics. The coex- 
istence of magnetic and electric order is associated to 
transition metal perovskites^ whereas spintronics pro- 
posals are based on materials where either spin-orbit 
interactions^ or d electrons^ play a prominent role. Here I 
show that the magnetic ground state predicted^'^i^ii'^i'^iiS 
in graphene zigzag ribbons, a chemically simple sys- 
tem without d electrons and negligible spin orbit cou- 
pling, is indeed a new electronic phase whose order 
parameter is the product of the spin and the elec- 
tric polarization. This new electronic state may be 
studied experimentally thanks to recent progress in the 
fabrication of graphen o^^d^d^ g^jjjj graphene based flat 
nanostrucutre a^^d^d^ . 

Early theory work predicts that graphene is a zero gap 
semiconductor, with electron-hole symmetry and linear 
conduction and valence bands. These features arise nat- 
urally from a tight-binding model with one Hz orbital per 
atom in a honeycomb lattice at half filling and they are 
related to fact that the honeycomb lattice is bipartite. 
The electronic structure of graphene nanoribbons de- 
pends dramatically on their atomic strucuturei^iiiii^ii^. 
Here I focus on graphene ribbons with zigzag edges. 
The single-particle description of this system features 
two almost degenerate quasi-flat bands at the Fermi 
energyi^iiiii^. These flat bands are associated to edge 
states . When Coulomb repulsion is added to this pic- 
ture within a mean field Hubbard model^i^i^ local mo- 
ments of oppposite signs form in the edges, with a total 
zero spin, and a gap opens at the Fermi energy. The 
predictions of this model are robust with respect to the 
addition of more orbitals, second neighbour hoppings and 
long range Coulomb interactions, all present in density 
functional (DFT) calculations^'^^-^-i". DFT results and 
the mean field Hubbard model yield very similar results 
for the low energy sector of the electronic structure both 
for zigzag graphene ribbons^" and nanoislands^-. This 
permits a significant computational simplification as well 
as conceptual advantage through the use of exact results 
valid for the Hubbard model^^. 



In this paper three things are done. First, I show that 
the mean field wave function of the Hubbard model for 
graphene ribbons with zigzag edges is that of two phase 
locked BCS condensates of spin-polarized electron-hole 
pairs living in the edge bands, the only bands affected by 
the interactions. Second, the BCS electron-hole coher- 
ence implicit in the wave funcion is associated to the ex- 
istence of spin-resolved transverse electric polarizations 
that yield a zero total electric dipole and spin when 
summed. Therefore, the standard mean field magnetic 
phase of zigzag graphene ribbons^!^'^'^'^'^ is an exci- 
tonic insulator phase with a hidden ferroelectric order. 
Third, the mean field bands are written in terms of the 
BCS gap and diagonal self-energies. 

Zigzag graphene ribbons are described with a single- 
orbital tight-binding mode l^^d^ plus a on-site Hubbard 
repulsion treated in the mean field approximation at half 
filling4,^>2i: 

H = ^ tfrpci^cf'a + U^np^-i{np^i) +nf^i{np^^) (1) 



where ct^ creates an electron at the H^ orbital of atom 

located at r with spin cr, np^a- = ct^cpa is the occupation 
operator. The first term in the Hamiltonian describes the 
first-neighbour hopping (t = 2.5eV) in the graphene rib- 
bon and the second describes on-site Coulomb repulsion. 
I take U ~ 2eV. The zigzag ribbon is a one dimensional 
crystal whose unit cell, shown in fig. la, is repeated along 
the X direction. The position r is determined by a unit 
cell index, x and a intra-cell index /. Notice that the top 
and botton atoms belong to different sub-lattices. In a 
unit cell there are with Np pairs of A and B atoms, and 
the width of the ribbon is ~ V3(iV)a, with N = 2Np. 

The spectrum of the self-consistent mean field Hamil- 
tonian for a ribbon with N atoms per unit cell has N 
bands per spin channel, half of which are occupied. Fig- 
ure lb shows the well known non-interacting (U — 0) 
bands for a ribbon with N = 18. Solid (dashed) lines 
represent full (empty) states. The two fiat bands in the 
outer region of the Brillouin zone correspond to states 
that are localized in the edges of the ribboni^. As shown 
in figure Id, they are not really degenerate except for 
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ka = ±7r. At zero temperature the lower (valence) band 
is full and the upper (conduction) band is empty. The 
operators that annihilate an electron in those bands are: 
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where W{c.v).k{I) ^-re the Bloch eigenstates and L 

is the length of the ribbon. The gap between these bands 
is proportional to the penetration of the edge states to- 
wards the bulk regioni^. 

Figure Ic shows the mean field interacting bands, 
shifted rigidly by — [//2. A gap opens at the Fermi en- 
ergy, in agreement with DFT calculationsi*^!^!^. The 
average spin-resolved charges along a unit cell are shown 
in figures 2a, 2c. Spin up (down) electrons pile at the top 
(bottom) edge of the ribbon and leave a charge deficit 
in the opposite side, also in agreement with DFT calcu- 
lations. Therefore, the edges have local magnetization 
with opposite sign. For a given spin, there is an excess of 
electrons in one edge that are missing in the other. The 
total electronic charge turns out to be the same in all the 
atoms. 
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FIG. 1: (Color online), (a) Zigzag ribbon unit cell, (b) Bands 
of iV = 18 zig-zag ribbon with (7 = (c) and U = 2eV. Only 
half of the Brillouin zone is shown, (d) Comparison of low 
energy bands. Interacting bands have been shift downwards 
by U /2. (e) U = Q single particle gap at the edge states. 



It is crucial to realize that the non-interacting and the 
shifted mean field bands are identical except for the low- 
est energy empty band and the highest energy occupied 
band which differ in the outer sector of the Brillouin zone, 
shown in figure Ic. These bands are denoted by — and + 
and V and c for the interacting and the non-interacting 
case. It turns out that both — and -I- can be expressed 
as linear combinations of c and v only, with an accu- 
racy better than 99%. Thus it is possible to relate the 
two interacting states — and -I- with the non-interacting 
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FIG. 2: (Color online), (a) and (c): spin revolved occupation 
nia as a function of vertical position in unit cell for N — 22 
ribbon, (b): Occupation factors Va{k)^, Ua{k)^ and ^^(fe)^ + 
Ua-{k)^. (d) Spin resolved dipole Vya{k) 



conduction and valence band through 
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Importantly, the spin dependence is limited to the 
phases. The moduli of the coefficients \uk\^ and \vk\^ , 
are shown in figure 2b. We find that vl + ul > 0.99. 
Using this relation we formulate the interacting theory 
in terms of electrons and holes in the non-interacting 
valence and conduction band. The mean field ground 
state, which is formed by filling all the mean-field bands 
below the gap, is written as |$) = |<I>)| x |$)| where: 
|$)ct = Ilfe/lj, ^|G"), where \G') denotes the state where 
all the bands below V (or — ) in figures lb and Ic are full. 
Making use of eq. ([3]) I write: 
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where |G)o is the non- interacting ground state with no 
holes in the valence band and no electrons in the conduc- 
tion band. Equation ^ is one of the important results of 
this work: the mean field ground state implicit in eq.p]) 
that yields the bands of fig Ic and the density profile 
of fig. 2c and 2d can be written as the product of two 
BCS condensates of spin polarized electron-hole pairs. 
This wave function is found in the context of excitonic 
insulator^SS, and non-equilibrium cxciton condensates^. 

Importantly, this BCS state implies the existence of 
non-magnetic long-range order for the interband oper- 
ators. The numerical calculations systematically show 
that the interband coherence 
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are finite and their relative phase is locked: i(l>l~4'l) — 
Interband coherence is zero for f7 = and is related to 
observables that mix the valence and conduction band. 
These acquire an anomalous expectation value in the 
?7 > phase. The interband coherence is associated 
to electric polarization^^ in non-equilibrium exciton con- 
densates and to electronic ferroelectricity in the case of 
Bose-condensation of slave bosons in the case of mixed- 
valence compounds^^. Hence, I look for the connection 
between the interband coherence implicit in eq. ([?]) and 
the spin-resolved electric dipole implicit in figs. 2a, 2c. 
The electric dipole is written as the sum of the spin- 
resolved dipole Vy — Vy^ + Vyi whcrc 

'Pyy ^'^yi(^'<^x,i,a = d^yGl^^^^Cky,a (6) 



xJ 



are the spin-resolved components of the dipole operator, 
and the dipole matrix elements are given by 



(7) 



The labels v, v' run over the non- interacting bands. The 
ribbon is centered at y = Q. Because of the mirror sym- 
metry of the zigzag unit cell, d^^(k) — 0, so that only the 
band-mixing terms on eq. jT]) can yield a contribution. 
The average dipole operator in the state ^ is: 

CPy^ = 7 E dcvik)u,vle-^*' + h.c. = jJ2 VyAkm 



Whereas {Vy) = the spin resoved components Vyaik), 
shown in fig. 2d are finite and with opposite sign. A zero 
net dipole resulting for the sum of two opposite spin- 
resolved dipoles is expected from inspection of figures 2a 
and 2c and the homogeneous spin-summed charge dis- 
tribution. Thus, the spin resolved dipoles are related to 
the interband coherence and the absence of net electric 
dipole is related their phase locking in eq.(l5])- In order to 
characterize this new kind of electronic order, I introduce 
the spin dipole operator: 



'{y,z) 
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where ^, is the Pauli matrix. Thus, the relevant or- 
der parameter associated to the electronic state (|4]) is 
Tra{^\'P{y, z)\^). Spin rotational invariance permits to 
choose z along any direction in the spin space. This order 
parameter is invariant under the combined action of time 
reversal and mirror symmetry, and provides a natural 
explanation to the spin-polarization of the system when 
subject to a transverse electric field, predicted by DFT 
calculations. Notice that this phase is different from the 
non-magnetic ferroelectric phase predicted IuIg!, which is 
not found in DFT. 

The mean field state invites to write the interacting 
bands in terms of a BCS-like gap related to interband 



coherence. To do that, I project out all the bands except 
C and V: 



(10) 

The occupation of the sites is expressed as nia = 5 + 
ami, where a — ±. This automatically ensures that the 
occupation in each site is 1. Using transformation eq. 
(fTU)l . the mean field Hamiltonian reads: 
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with 



U 



^ua{k) = e,(fc) + - + C/a^|7/-fc,,(/)|2(m,) (12) 



where e„{k) are the U — bands and v = c,v. The 
second term in (IT2]) is the rigid shift of the bands 
and the third term is the diagonal self energy Y,^ cr{k). 
The off-diagonal self-energy reads: 



mi) 



(13) 



Notice that Acr(fc) = ~Aa-{k), which explains the phase 
locking of eq. ([5]). Notice also that in the Hubbard 
model the self energies for spin a electrons depend on 
the density of carriers with opposite spin a. For each 
ka the mean field two by two matrix can be written as 
+ hcr{k)f where r are the Pauli matrices, and the ef- 
fective field can be written as: 

Kik)= (Re{A^ik)),ImiA^{k)),^^^^^^^^) (14) 



The eigenvalues of this two by two matrix are 



1 



E±Mk) ^2^^^ V (?-.- - ^vAkW + 4|A,(fc)|2 



The transformation ^ permits to diagonalize (jlip . ob- 
taining n = J2k,cr,r=± Er,a{k)jl^„fkTa- At zcro temper- 
ature only the lower branches E-^„{k) are occupied. The 
mean field dispersion Era{k) depends on ^^aik) and on 
A^{k) which in turn depend on the magnetization: 



m{I) ^J2^kAl)i'k,v{I)ukvl+h.. 



(15) 



The magnetization depends on the transformation fac- 
tors, u and V, which depend on the energies through: 



\vk\' = lU- 
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Equations (|12I13I14I15I16P form a self-consistent set. 
The numerical solutions of the mean field Hubbard model 
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FIG. 3: (Color online). Relation between dispersion E+cr{k) 
and the diagonal Eyo-(fc) and off-diagonal A(fc) self-energies 
for two ribbons with N = 22 (a) and iV = 42 (b). 

(P) also satisfy these equations. This permits to re- 
late the mean field dispersion to the diagonal self energy 
\'^c,a{k) — T,y^cr{k) \ and the off diagonal self energy |A(fc)|, 
shown in fig. 3. They are spin indcpedent. It is apparent 
that these self-energies are finite in different regions of the 
Brilloiun zone. The diagonal self-energy are related to 
the non hybridized edge states located in the outer region 
of the Brillouin zone whereas the off-diagonal self-energy 
occurs for weakly hybridized edge states, at smaller 
Further reduction of |fc| opens the single-particle gap, 
which overshades the self-energies. These results also 
permit to unveil the origin of the gaps and A'' in- 
troduced in d (see fig. 3). It is apparent that the A^ 
gap is given by the diagonal self-energy whereas the A" 
gap is related to the non-diagonal self-energy. Accord- 
ingly, A^ is insensitive to the ribbon width and can be 



approximated by A-'^ ~ 2?7|m| where |m| is the magne- 
tization of the edge atoms. In contrast. A" decreases as 
the ribbon width increases due to the smaller inter-edge 
hybridization. 

The long range order implicit in this and previous 
mean field theories of graphene zigzag ribbons^i^iiiSiiSiiiS is 
known to be destroyed in one dimension because of long- 
wavelength spin wave modes^ associated to the breaking 
of a continous symmetry. I have verified that the results 
of this work remain valid for finite length graphene rib- 
bons and tubes for which Goldstonc modes have a con- 
finement gap. Therefore, the results of infinite ribbon 
systems are relevant for finite systems. 

In summary, the BCS wave function ^ describing two 
phase-locked condesates of spin-polarized electron hole 
pairs is the the collective wave function behind the in- 
sulating ferrimagnetic phase in graphene zigzag ribbons 
portrayed by mean field Hubbard'*''^'^ model, which yields 
the same results than DFT calculations- The un- 
derlying electron-hole coherence in each spin-channel is 
related to mirror symmetry breaking of the charge den- 
sity for a given spin. Their relative phase-locking war- 
rants that the total spin and electric dipole are zero. 
The natural order parameter for this electronic state with 
magnetic and spin-hidden electric order is the spin-dipole 
operator (eq. ([9]). The reformulation of the mean field 
theory in terms of a 2-band BCS model rationalizes the 
shape of the mean field bands in terms of diagonal and 
non-diagonal self-energies. The joint presence of electric 
and magnetic order anticipates non-trivial magnetoelec- 
tric effects in graphene zigzag ribbons. 
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